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I. INTRODUCTION 


Understanding of the characteristics of the airflow over 
an airfoil is of paramount importance CoO.esEne Vasntoil) 
designer. Two methods are currently available which give 
accurate results. The first is the use of wind tunnel tests. 
The drawbacks to this method are cost and time consumption. 
The second is the processing of the Navier-Stokes equations. 
This method's drawbacks are the requirements and expense of 
using supercomputers due to the extensive calculations and 
storage requirements. There is still a need to come up with 
an inexpensive, fast and accurate engineering tool to compute 
airfoil flows. 

Several methods have been derived to accomplish this end. 
But the most promising is the Viscous-Inviscid Interaction 
method. The outer flow is computed using inviscid flow 
equations, and the inner flow is computed using Prandtl's 
boundary layer equations. The key to this method is the 
extent of interaction between the inner and outer flows. 

The purpose of this thesis is to evaluate the capability 
of the viscous-inviscid interactive aircode developed by 
Tuncer Cebeci and associates at the Douglas Aircraft Company 
[Ref. 1]. This computer program was applied to four airfoils 
with various angles of attack and Reynolds numbers. The 
computer results were then compared to previously reported 


experimental results. 


The conservation of mass and momentum are summarized in 
Chapter 2, inviscid flow calculations are discussed in 
Chapter 3, and viscous flow equations are described in 
Chapter 4. Viscous calculations are presented in Chapter 5, 
and the specific interaction methods are shown in Chapter 6. 
Finally, in Chapter 7 computer and experimental results are 
compared for the NACA 663-018, 0010 (Modified) and 4412 


airfoils as well as the Wortmann FX 63-137 airfoil. 


II. FUNDAMENTAL EQUATIONS 


The conservation of mass and conservation of momentum 
provide the foundation for incompressible flow analysis. 
With these fundamental concepts along with appropriate 
assumptions and approximations, working relations for two- 


dimensional, incompressible flow are obtained. 


A. CONSERVATION OF MASS (CONTINUITY ) 

The conservation of mass principle states that mass 
cannot be created nor destroyed. Equating this statement to 
a fixed control volume the net mass flow rate into and out of 
the control volume equals the time rate of change of mass 
within the control volume [Ref. 2:p. A-1]. 

Given a control volume, the mass flow rate through one of 
its surfaces is equal to the product of the fluid density, 
the fluid velocity normal to the surface and the area of the 
surface [Ref. 3:p. 29]. 

dad mass 
———— = Vens L2aa) 
dt 
In 2-D flow the x-component of the mass flow rate at the 
center of the positive x-face, position dx/2 and side length 


dy, is represented by Taylor series expansion [Ref. 2:p. A-2] 


O dx (¢? ax 1 
pu + —(pu)— + Cg tt. Ye (2 2p 
Ox 2 Ox? 2 a 





As ax approaches zero all higher order terms disappear 


leaving 


ou + —(pu)— dy. (27235) 
Xx 2 


Similarly, the x-component of the mass flow rate at the 
center of the negative x-face, position -dx/2 and side length 


dy, is represented by 


©) dx 
pu te ars dy. (2.4) 
x Z 


As illustrated in Figure 2.1 the net mass flow through the 


four sides of the 2-D control surface is 


a) dx 9) ax 
pu - —(pu)— jdy -{] pu + —(pu)— [dy 
Ox 2 Ox 2 
dy dy 
+) pv - 9 ov) dx -| pv + ly gy ax (2.59) 
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Figure 2.1 Mass Flow Through 2-D Control Surface 
[Rete 4:p. 12} 


which is equal to the rate of change of the mass within the 
control volume 


OP 
= OXY 5) 
O 


t 


Combining (2.5) with (2.6) and simplifying yields 


- 2-( pu) axay - OF mages = oP axay (Zama) 


x OY x 


Dividing by dxdy and rearranging yields 


OP + Ci + OF run = 0 (228) 


Or Ox OY 
For steady, incompressible flows the continuity equation 
becomes 
Ou Ov 


ne gr 0 ( 23335) 
Ox Oy 


Or in vector form the continuity equation [Ref. 3:p. 30] is 
Vee Ve=0 (220) 


B. CONSERVATION OF MOMENTUM (NAVIER-STOKES ) 
The conservation of momentum, Newton's second law of 
motion, states that the rate of change of the linear momentum 


is equal to the sum of the forces applied [Ref. 2:p. B-1]. 


Tee ( mV) (el 


As illustrated in Figure 2.2 the two significant forces which 
act on an element of fluid are surface forces which act on 
the surface only, pressure and shear, and body forces which 
affect the mass of the element, such as gravity. Assuming 
moment equilibrium in an element, Txy = Tyx, the 2-D first 
order Taylor series expansion for normal and shear surface 


forces in the x-direction is 


dx a dx 
Vsene + = Ceo) -— Uxx + — (Tx) — dy 
x 2 Ox 2 


dy 4) dy 
Uxy + E(t) — -— Txy + — (Tx y)— dx 
2 OY 2 


= Om evaay + Oe. \axay: (2 tce} 
Ox OY 


The body forces per unit mass are represented by 
Fsopvy = Xi + Yj + Zk (22a) 
such that the x-component of the body force on an element is 
fxcpony) = pdxdy*1X. (2.14) 


Combining equations (2.12) and (2.14) the sum of the forces 


in the x-direction is 


dy 
Tyy + (ae 


OY 2 
dy 
OY 2 
a 
ax 
a dx 4 txy + —( 2 
a 
dy x 
| ax 
2) dx —S  Txx t—( Tx) — 
toa ax Ox 2 
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Figure 2.2 Stresses on a 2-D Control Surface 
[Ref. 4:p. 15] 


SFs =| PX + ——(Txxn) + —(Txy) | dxdy. 2.15 
Ox Oy . | 


The rate of change of the linear momentum in the x-direction 


assuming constant mass is mdu/dat. 


du Ou dx Ou dy Ou Ox 
ee = —— — + — — + — Yila the chain rule, and — = u, 
at Ox dt Oy dt Ot | Ot 
OY 
and — =v the x-direction change in linear momentum for 
t 


particle is 


du Ou Ou Ou 
m— = pdxdy(v— + v— + —). (2216) 
dt Ox Oy oO 


Substituting equations (2.15) and (2.16) into the x-component 


of equation (2.11) yields 


Ou Ou Ou | D 
Gemeay {U—— + v— + ——) = Cant (Cao leet ——( Cs.45)  |AXdyY. (2.17) 
x Oy Ot x Y 


Now, in order to have the entire equation as a function 
of velocity the normal and shear stresses must’ be found in 
terms of velocity. 

By assuming a Newtonian fluid [Ref. 2:p. B-5] the shear 
stress is linearly related to the rate of angular deformation 


with fluid viscosity being the proportionality constant. 


‘0 


t= (23s. 


The normal stresses are equal, but opposite in direction 
to the pressure when no shear stresses are involved. With 
Shear stress from viscosity it is assumed that the normal 
stresses deviate from =P and that the deviation is 
proportional to both a) the rate of linear strain in the 
Girection concerned, andb) the rate of volume deformation. 
Therefore, the normal stress in the x-direction [Ref. 1l:p. B- 
10] is 

Ou 2 OO 
Uses = Rt Ci a ee | ae 
Ox 3 Ox Oy 
Applying the conservation of mass, equation (2.9), equation 


(2.19) simplifies to 


Tax =O ome beck — (2.20) 


Substituting equations (2.18) and (2.20) into (2.17) and 


dividing by dxdy yields 


Ou Ou Ou 3) Ou Oo ou Ov 
(U— «+ V— + —) = OX + —(=-P + 2n(—)) + D—(— + =i 
j a Oy oat i Ox : Ox = Oy Ox 


After multiplying and rearranging the right hand side becomes 
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OP O7u O7u O2v 
PX - —— + Qu ap + 
Ox Ox? Oy? OYOX 











which is also equal to 


OP O2u O2u Oo Ou Ov 
-—— tH +H + p—( 
Ox Ox2 Oy? Ox Ox Oy 


——— 


px 








Again applying the conservation of mass, equation (2.9), 


equation (2.21) becomes 








Ou -Ou Ou OP O2u O7u 
pCu—— + V—— + —) = OX = — + TI + - 2tee2 
Ox Bye Ole Ox Ox2 Oy? 
With v= yp/p = kinematic viscosity and neglecting the body 


force, X, the two dimensional Navier-Stokes, conservation of 
momentum equation for incompressible and constant viscosity 
flow in the x-direction [Ref. 2:p. 436] is 

Ou Ou Ou Orr OZu Ou 


8 eS ee a | + 


Ot Ox Oy 9 Ox Ox2 dy? 








). Gro ) 


Similarly, in the y-direction the Navier-Stokes equation is 


Ov Ov Ov 1 OP O2v O2v 
— + U— + V— = - — — + + 


ot Ox OY 0 OY Ou By- 


ye (2.24) 
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IIIT. STEADY INVISCID FLOW 


Although real fluids are viscous the major effects of 
viscosity are concentrated in a region, or layer adjacent to 
a body. Therefore, analyses of inviscid flow are useful and 
serve as a good approximation to flow outside the boundary 
layer and wake behind the body. 

The justification for applying the results of perfect 
fluid analyses to viscous flows was postulated by Ludwig 
Prandtl in 1904 [Ref. 3:p. 299}. He stated that the effects 
of viscosity on the flow around streamlined bodies at high 
Reynolds numbers are effectively limited to a "thin" boundary 
layer. The characteristic length to judge thinness is the 
distance from the forward stagnation point to the point of 


consideration. 


A. VELOCITY POTENTIAL 

For flow outside the boundary layer it 1s a great 
advantage to Simplify equations and develop aesingle 
governing equation. With the assumptions of steady flow, no 
energy transfer to or from the fluid, no body forces, no 
shear stress (inviscid), and irrotational flow the velocity 
potential, 6, 1s wtilizged [Rett 3p. ecole 6, a scalar 
function of spatial coordinates, x and y, is defined such 


that 


V = Vo (3 a8 


eZ 


and 
i= == Vv=— B--2) 


The importance of the velocity potential is that only one 
equation is needed to describe the irrotational flow. 


Velocity components are obtained using equation (3.2). 


B. LAPLACE EQUATION 

For steady, incompressible flows the continuity equation 
(9) is 

Ou Ov 
SS. + eae 
Ox Oy 
Rewriting (2.9) in terms of the velocity potential the 
equation becomes 
O°> 07 


+ = O (e323 ) 
Ox, Ove 








This form of the Laplace equation [Ref. 2:p. 81] is the 
governing equation for steady, irrotational flow of an 
incompressible fluid. 

The importance of equation (3.3) is that it is linear 
allowing for the principle of superposition. For example if 
$., dO2, ds...are solutions of (3.3), then the sum @¢@ = $1 + do 
+ d3 +...is also a solution of (3.3). Superposition of 
lrrotational, incompressible flows allow for the construction 


LS 


of complex flows that are also 1rrotationa: and 


incompressible. 


C. SIMULATION (CONFORMAL MAPPING) 

The inviscid flow about an airfoil can be obtained most 
conveniently by means of a transformation, which starts with 
a known flow about a simple contour, a circle, distorts the 
contour into the desired shape, and simultaneously adapts the 
flow to that shape. The transformation is accomplished using 
a sequence of three conformal mappings [Ref. 1:p. 324]. 

The first mapping, necessary only when the airfoil 
trailing edge has non-zero thickness, is accomplished using a 
logarithmic mapping function. The airfoil is perturbed 
slightly to make the upper and lower surface, trailing-edge 
points coincide. 

The second mapping analytically removes the trailing-edge 
corner using the Karman-Trefftz mapping. 

The third and final mapping transforms a quasi-circular 
shape into a perfect circle using an iterated sequence of 
Fast-Fourier Transform applications. 

During the transformation of streamlines about a circle 
to those about an airfoil, the preferable approach insuring a 
transformed flow free from vorticity uses the complex 
Variable z = x +iy [Ref. 5:p. 285]. The transformation of z 
to another plane is E= f(z) = & + in. The potential 
function, Q = 6+ i¢, is irrotational and incompressible in 
both planes. 
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The streamlines, ¢, and equipotential lines, $, of a flow 
in the zZz-plane will transform into another orthogonal network 
of lines in the (-plane. Different magnification ratios and 
different amounts of rotation at different points in the 
field will, however, change the appearance of the flow 
pattern from about a circle to about the airfoil. 

The general transformation function whose derivative 
at/dz satisfies the requirement d{/dz approaches unity for 


large values of zis 


ae tere: tea /Z2 ENC. /22°455. | (3.4) 


The requirement is necessary as streamlines are not distorted 
a great distance from the body where the body's shape has no 
influence on the flow. 

The coefficients may be real, imaginary or complex. A 
finite number of the coefficients are determined from the 
specified normal velocity components equally spaced around 
the unit circle, and from the Kutta condition which ensures 
stagnation at the trailing edge. 

While in the first iteration the normal velocities are 
zero, and the solution for flow over a circle is used, the 
subsequent normal velocity boundary conditions are determined 
from the previous viscous-flow calculations using the 


equation 


its) 


d 
Vn = —(Uu.6”™) (3258 
ds 
where ue is the velocity at the edge of the boundary layer 
and 6* is the displacement thickness. Once the coefficients 


are found, the real and imaginary parts of equation (3.4) are 


equated yielding 
&= &(x,y) and N= (x,y). 


As x7 + y? = r? the two equations of & and n are transformed 


to 
xX = x(&,r*) and y = y(n,r*). 


Then x? and y? are added to yield 


x= + y* = 7r? = | (eer?) |? a | vine?) |? (3 768) 


After dividing both sides by r? 


al 1 
2 =— |} X(€,r7) |? + — | y(n,r*) |. (3.7) 


Then each circle of radius, r, in the z-plane is transformed 

to the proper shape in the (-plane to describe inviscid flow. 
1. Transformation of Velocities [Ref. 5:p. 291] 

In the z=-plane as Q(z) = o(x,y) + i¥{x,y) See 


velocities are defined by 
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adQ(z) 





= Wise ah 1b Ae eS ou 


where, 


In the (-plane the velocities are also found 


aQ 
— wie Vy eed V 3.98) 
at E 1 
where, 
Vie = Vg + 1Vy . 


The velocities in the two planes are equated by 


dQ az eee Ve 
VE - 1V =- —_—_—- — .) -lCIrh eeoesseererowmwrw"""_—ss (5G) 
az ac ac /daz 


The pressure in the transformed stream is related to the 


stream velocity through Kelvin's equation 


1 
2 + P = Constant = —pop 
2 


Z 
5 = 3.1.1.) 


V 
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Ue 
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IV. VISCOUS FLOW 


A. DERIVATION OF BOUNDARY LAYER EQUATIONS 

The previous analyses provide a valid solution to the 
flow outside the boundary layer. Within the boundary layer 
however, the effects of viscosity cannot be neglected. In 
laminar flow governing equations can be obtained by 
Simplifying the conservation equations. In turbulent flow, 
however, the number of variables outnumbers' the equations. 
Great dependence is then placed on dimensional reasoning and 
on hypotheses suggested by experimental results. 

The most important deduction from Prandtl's thin boundary 
layer theory is that static pressure can be considered 
constant across the boundary layer [Ref. 3:p. 299]. 


OP 
een — Tae @ (4.1) 


OY 
As the boundary layer thickness, 5, is small, d6&/dx is also 
small. Streamlines are then only slightly curved and the 
radii of curvature, R, are large. With a large R the 


equilibrium condition 


18 


ii@wetrates that oOp/ Oy will “be very small and can be 
neglected. Experimental results confirm that Op/ Oy may be 
neglected even over surfaces of small radii of curvature. 
Also, in a thin boundary layer with a slowly changing 
thickness Ov/Ox is much smaller than Ou/Oy. The significance 
is then that the normal shear stress may be neglected when 
compared with the viscous’ shearing stress. Equation (2.18) 


then becomes 


Ou 
Txy = H—. (4.2) 
OY 
With this simplification the approximate equation for 
flow in the two-dimensional boundary layer can be found 


directly. Newton's second law of motion applied to a fluid 


element of mass may be written 





Ou Ou Ou a} ro eee 
Oe ey) = (= j)dxdy (4.3) 
Ot Ox OY Ox Oy 
as illustrated in Figure 4.1. Substituting equation (4.2) 


and dividing both sides by dxdy yields 


Ou Ou Ou opm O eu 
p(—— + u—— + v—) 


( 
Ot Ox OY Ox oy OY 


i 
a 
ij 
n> 
a 


In terms of kinematic viscosity equation (4.4) becomes 


Ou Ou Ou 1 0p O7u 
= + (4.5) 
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Figure 4.1 Forces Acting on an Element in the Boundary Layer 
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This equation is the boundary layer equation of motion and is 
identical to the equation found using an order-of-magnitude 
analysis [Ref. 3:p. 443]. Equation (4.5) is also nearly 
identical to the Navier-Stokes equation (2.23) with the 
exception that the term y0O7u/Ox? is deleted. The order-of- 
Magnitude analysis suggests that this term, y0O7u/Ox?, may be 
neglected compared to yO2u/ Oy?2. Combined with the 
continuity equation Ou/Ox + Ov/Oy = 0 (2.9), equations (4.5) 
and (4.1) are known as Prandtl's boundary layer equations 
[Ref. 2:p. P-9]. 

For an incompressible flow, there are three variables, u, 
v and p, but only two equations, (4.5) and (2.9). The 
equations may be solved though, by first determining p as a 
function of x uSing inviscid methods, setting Op/Oy =) Oaien 
the boundary layer, and then solving (4.5) and (2.9) for the 


velocity distributions. 


B. TURBULENT FLOW 

Turbulent flow as differentiated from laminar flow is 
characterized by fluctuating instantaneous properties which 
greatly increase the complexity of the problem. A very 
useful simplification to the turbulent problem is then the 
use of time-averaged values, denoted by a bar over the value. 
Instantaneous values are indicated by the prime [Ref. 4:p. 


Zone 


Pa! 


The continuity equation containing total values becomes 


Cie one 
ma Ue 
Ox OY 


Simplifying the equation becomes 


oa _ ) Os) (v') 0 
—(u) + —(u') + —(v) + —(v') = 
Ox Ox OY y 
with 
oO _ = 9) O 
—(u) = —(u) and —(u') = —(u') = 0 
x x x Ox 
ee a qd (Va) (v' ) 0 
(Py an —(v') = —(v') = 0. 
Ox Ox Ox Ox 


The time-averaged continuity equation for turbulent flow is 


now 


Cs eee (4.6) 
Ox Oy | 


Ze 


Applying total values, the steady version of the Navier 


Stokes equation (2.23) becomes 


Ou + uy O- (uae uu") 
», 9 a ns i, (4.7) 
Ox? Oy? 


After simplifying, the time-averaged Navier-Stokes equation 


for turbulent flow [Ref. 1:p. C-10] becomes 

















Ot Dow) 1 OP Onl» Oem 

u-—(u) + Vv eee ny | + ) 

Ox OY 0 Ox Ox2 Oy? 

= (oe) =) -—(u'v’). (4.8) 
Ox OY 


The new terms, 0/Ox(u'2) and O/Oy(u'v'), which correspond to 
normal and shear stress, are called Reynolds 


stresses. The similar y-component terms’ are )/dy(v'2) and 


gu’ ) /Ox. 


C. TURBULENCE MODELS 

The time-averaged Navier-Stokes equation is nearly 
identical to the original equation except that the 
instantaneous values are replaced by the mean or time- 


25 


averaged values and two additional terms involving 
fluctuating velocities, u' and v', appear. An interpretation 
of these two terms compares them to the previously existing 
terms O7u/Ox2 and O7U/Oyz. The right hand side of equation 
(4.8) less the pressure term, and after multiplying by 


density, becomes 




















O-u oeu O 9) 
y + ea —— tj ' < nae _—_.¥y ! t 
Ox? "Dy? "Dy: i _ 
Or 
Ounou . Ou 
| re 12 ae | eu 
5x ae OE eas ee pu'v') 


As each term has the dimensions of stress, and p(Ou/Oy) is 
part of the laminar shear stress UTyx, it appears that the 
term -pu'v' represents a turbulent addition to shear stress 
f[Ref. 2:p. T-2]. Now, this shear stress is really a vertical 
mixing of horizontally, travelling fluid particles. A model 
of this mixing then calculates the rate of momentum transfer 
involved. 
1. Prandtl's Mixing-Length Model 

To predict the turbulent stresses Prandtl assumed 
that turbulent fluctuations are primarily the result of the 
mean velocity differences between two layers in the flow. 
Therefore, u' is proportional to Ou/Oy with the 
proportionality factor having the unit of length. 
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Ou 
ea) a — (4.9) 
OY 





Also, assuming that v' is of the same order of magnitude as 


u' at a particular point, 


Ou 
(v'2)# = pb. (4.10) 
OY 





Substituting u' and v' into the turbulent shear stress, tT-, 


is 


Ou Ou 
- opu'v' = - pgab(—)(—). (4.11) 
Oy 0 





As a andb are both unknown constants of length, they both 
may be replaced by the "mixing length", 1, the hypothetical 
distance between the two layers involved. The turbulent 


shear stress is now 


Ou| Ou 
| a (4.12) 
Oy| Oy 


tyr = ol? 





which insures the correct sign. 
2. Eddy Viscosity - Cebeci-Smith 
The turbulent addition to shear stress may also be 


modeled in terms of "eddy viscosity". As laminar shear 


ao 


stress = lyx = n(Ou/Oy) = pv(Ou/ Oy, turbulent shear stress 


may be equated 





Ou 
“= pu'v' = PVe—, (4.13) 
Oy 
where yet, the eddy viscosity is empirically determined. 


The method used here is that of Cebeci and Smith as 


modified by Cebeci, Clark, Chang, Halsey and Lee, [Ref. Wee 














BZ ic 
Eddy viscosity by Cebeci is represented by 
4 u 
{0.4y[1 - exp (- —)]}? |—J Yer (0 S$ Y S Ye) (4.14) 
A Y 
the = 
a | (Ue - u)dy | Yer (Ya < Vyas 30) (4.15) 
oO 
where 
26 i 
A= ——— and a 
3 fot 5. Siliv 76) = 
u 
Yo 
Y 


The distance from the body y. which is less than the boundary 
layer thickness, 6, is the distance where the two equations 
(4.14 and (4.15 give the same resultant V:-. 


26 


The intermittency factor, Yer, Which indicates the 


local fraction of turbulent flow to total flow, is given by 


me ax 
Yer = 1 =- eCxXp [= G(X -—"X-L) a a (4.16) 
Rer Ue 


The location of the start of transition iS Xexr, and the 
empirical factor G is 
1 ie 


CO R- ame (A 7) 
120027 ae 





where Rxtxr is the transition Reynolds number 


Ree = (Wax /V es. (4.18) 


In equation (4.17) the empirical constant GYer = 
1200, was used by Chen and Thyson [Ref. 4:p. 327]. Values 
lower than 1200 may give better results at low Reynolds 
numbers as will be discussed in Section VII. 


The expression for a is 





0.0168 0.0168 
i (4.19) 
| 205 
Ou/Ox 
1-8 
Ou/Oy 


The non-dimensional factor F represents the ratio of 


the product of the turbulent energy by normal stresses to 


EG 


that by shear stress evaluated at the location where shear 


stress is maximum. 








(u'2 - v'2)Qu/Ox 
EF Se) = oe (4.20) 





- u'v' Ou/dy 


=. —(- u'v') max 





According to the data of Nakayama [Ref. l:p. 327], 8& can be 


represented by 














6 
ees |) FC Cn 
Mo saree 1 + 2Rr(2 - Re) 
R= oe = 
as Ley 1 a5 R- 
l= U'V' ) wemx ——$—— Rr 2 1 
Rr 
where Rr = Tw/(- wu'v')maxw and tw is the wall (body) shear 
Stee coe 


D. TRANSITION 

The location of the onset of laminar-to-turbulent 
transition when not found experimentally is determined 
empirically. The method used by Cebeci [Ref. 4:p. 333] is 
the criterion proposed by Michel. 

At the point of transition the Reynolds number based on 
momentum thickness, ©, is related to the Reynolds number 


based on the coordinate position, x. 
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Retr 


iverece x aS ateZ} 


Rextxr 


where 


Rea = Ua(X/y ) and Retr = U.(O/y). 


In the Cebeci Code the transition may be determined in 
the following ways: 


1) The points of transition are calculated using Michel's 
ei) tel ron . 


2) If laminar separation occurs forward of the criterion 
points, Michel's criterion is disregarded and 
transition is redefined at the separation point. 


3) The transition locations may be specified by the user. 
provided stall is not computed. 
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V. VISCOUS METHODS 


Momentum transfer in fluids is accomplished by 
hydrostatic pressure and viscous’ stresses. When viscous 
stresses are negligible, fluid behavior can be predicted by 
inviscid flow methods as stated in Section III. 

Viscous stresses caused by a variation in velocity ina 
direction normal to the flow are called shear. The most 
common shear is that found in the boundary layer between a 
displayed stream andthe solid surface. With the "no slip" 
condition fluid velocity is zero on the surface, but the 
velocity gradient is not so constrained. From the body along 
its normal direction the fluid velocity asymptotically 
approaches that of the free stream. 

As mentioned in Section III, Prandtl hypothesized the 
division of the flowfield into the two regions, the boundary 
layer where viscous effects cannot be neglected and the 
region outside the boundary layer which may be considered 
inva Ged: 

This hypothesis allows for the use of the parabolic 
boundary layer equations of section III instead of the 
elliptic Navier-Stokes equations. Depending on the boundary 
conditions, solutions fall into three methods [Ref. 6:p. 13]: 


1) The direct boundary layer method. This method uses the 


“nO” Skip" condition, where normal and tangential 
velocities are zero at the surface, and a 
pre-determined velocity 1S specified at the boundary 
layer edge. 
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2) The inverse boundary layer method. Boundary conditions 
are replaced by wall shear or displacement thicknesses. 


3) The interactive boundary layer method. The edge 
boundary condition drives a combination of displacement 
thickness and external velocity. 

Methods one and three will be discussed as they apply to 


the Cebeci Code. 


A. DIRECT BOUNDARY LAYER METHOD [Ref. 6:p. 13] 

This method for solving the boundary layer equations is 
used only near the leading edge where the viscous effects are 
small. Initial conditions are generated at the stagnation 
point and the equations are integrated around the leading 
edge. The numerical solution utilizes a finite difference 
method where the continuity and momentum equations’ are 
redefined as a system of linear algebraic equations. 

The method begins by describing steady, incompressible, 
2-D flows in a curvilinear coordinate system where x is 
directed along the airfoil surface and y is perpendicular to 
x. The boundary layer equations with the turbulent Reynolds 


stress are 





Ou Ov 
sa = 0 (oe k) 
Ox saOy 
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2 haa = + ay 5.2) 
prox.  pmOy. ae i = 5a 


a 


where the order of magnitude of any turbulent stress is 
assumed to be that of its laminar stress. The boundary 


conditions are: 
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The eddy viscous stress previously defined is reprinted as 





Ou 
-~pu'v!’ = ao (4 Gily 
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Also, the pressure gradient term may be written 


1 Op Oue 
a 2S ue—- (5.6) 
0 Ox Ox 


Therefore, the momentum equation (5.2) may be rewritten as 


du +00 due Od Ou 
u— V———— = eee ee | (De) (5279 
Ox Oy Ox Ox OY 

where b = 1 + yve/y and the boundary conditions are 


at ye= O° Wis. 0) f=on vseano) 20 (5.8) 
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1. Falkner-Skan Transformation [Ref. 6:p. 14] 

To solve the new boundary layer equations the 
Falkner-Skan transformations are used, which reduce the 
number of variables, and scale the normal coordinate y and 
the stream function $ with reference to the external 


velocity. 


Ue | 3 
n=y}— one) 
Vx 
Pale Uae EET) (5 1G) 


The continuity equation is automatically satisfied 
using the stream function (u = O¢/Oy and v = -oO¢/Ox). 
Therefore, only the momentum equation needs to be solved, 


which after transformation becomes 





m+ 1 OF ' Of 
fot") a+ fio leone | = x( fio = FN )(5.11) 
2 Ox Ox 
where m = (xX/u.)(due/Ox), a dimensionless pressure-gradient, 
and f' = Of/On. 
This equation (5.11) is a third order partial 


differential equation, and the solutions are "non-similar" as 
they are functions of both x and yn. If the solutions were 
mye a function of 7, then the right hand side of the 


equation would equal zero, andthe flow would be "similar" 
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[Ref. 1: Dp lav—1Colr To solve equation (5.11), a numerieee 
solution is needed such as the "box" method. 
2. The Box Method 

The box method, developed by Keller in 1970 [Ref. 
1:p. 331), 1s a widely used methods for solving non-linear 
differential equations. The steps of this method ineimee the 
conversion of the Falkner-Skan, transformed, momentum 
equation into a system of first-order partial differential 
equations. This non-linear system, after conversion from a 
continuous function to discrete, is linearized by Newton's 
method. The block elimination method is then used to solve 
the linearized difference equations of the boundary layer 
problem. 

The third order momentum equation (5.11) is converted 
into a first order system with the addition of the dependent 


variables U and V [Ref. 6:p. 14]. 





U = £! (5. i 

v= vU' = £" ( Seiese 
m+ 1 OU Of 

(bV)' + £V 4+. m(ie = Us) 233 (SV) (5.14) 
2 Ox Ox 


The boundary conditions are 


at 7 O; Ux, 0) = 0; brieo) = 6 


at N=Ne; U(xX,Ne) = 1 
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The solution domain, O<xs<xxz and O<n<ns, of the 
continuous functions f£, U and V is then covered by a 
rectangular grid to facilitate the problem solving with a set 
of discrete values. This grid is shown in Figure 5.1. 

The subsequent notation [ ]4 is used to represent the 
quantities of £f, U or V at the point (x:, n;3). 

All quantities can then be approximated by network 
coordinate values. Using the grid system, the solution of 
the parabolic layer equation at a certain streamline position 
depends solely on the solution of upstream positions, while 
no downstream influence needs to be considered. As 
calculations proceed from the stagnation point in the 
downstream direction, the overall solution can be obtained 
incrementally. Hence, one step of the solution procedure 
sets up the governing equations for a column of grid boxes in 


the sub-domain 


Xi-18XSXi and OSNSHG 


and solves for the values of the downstream grid position. 
The x-grid position currently solved for is then assigned the 
superscript "i" while "i-1" represents the previous position 
of known properties. Using coordinates of box midpoints and 
centered-difference derivatives, the equations are actually 
satisfied midway between the grids. 

Equations (5.12) through (5.14) in terms of finite 


differences [Ref. 6:p. 15] are now written 
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Figure 5.1 Rectangular Grid for Finite Difference 


Approximation 
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where the ordinary differential equations, (5.15) and (5.16), 
are centered about (xi, ns-%) and the partial differential 
equation, (5.17), is centered about (X:si-34, Ns-a). 

The boundary conditions are 


= gare oe 
= Oy fi = 0 and US = Jl 


Equations (5.15) and (5.16) are the centered difference 
derivatives. 
3. Newton's Method [Ref. 6:p. 15] 
This set of finite difference equations is nonlinear 


with combinations of unknowns. Newton's method is therefore 


Sy 


needed to solve the system. The variables are linearized 


using an iterative procedure with preceding values. 
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K is the iteration counter. After substituting these values 
into equations (5.15) through (5.17) and neglecting the 

1,K Kk 1.K 
quadratic terms of 6f OU and 6V , the system of 

a 4 J 


unknowns is then linear. 
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The boundary conditions are 


5fo = OF 65Uo = 0, 5U 5 = 0, 
4. Block Elimination Method [Ref. 6:P. 17] 
The system of equations are iteratively solved until 

1K i.K 1.K 
6f ,2ouU and 6V become small enough to be neglected. 

5 7 : 
The solution method is that by Keller and is called the Block 
Elimination Method. In this method block-tridiagonal 


matrices are composed of blocks. Only those blocks on the 


main diagonal and on the two adjacent diagonals are non-zero. 
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where the blocks are 3x3 matrices 
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The right-hand sides of the equations are 


1k 


i.K-1l i.K-l 
ite i ae = “tS + h.uU. 2 Shi 
( 1) 5 j-1 j j°j-3 
Bays" - {As defined in equation (5.20)} 2<5<I 
a ee 5 1.K-1 i.K-1 1.K-1 
zi z <h< 
(ra) U, Us iy + Ne 543 ech su 
alae Soe ijk digg wee 
(Tr), a OF (C2), on OF (Ts) 5 = 0 
The block elimination method solves the linear 
equations with two steps. The forward step eliminates the 


lower diagonal of the tridiagonal matrix. The reverse step 


solves the remaining system from bottom to top. 


B. INTERACTIVE BOUNDARY LAYER METHOD [Ref. 6:p. 18] 

As the direct boundary layer method 1S, as previously 
stated, restricted to regions of small viscous effects, and 
integration of the boundary layer equations fails at points 
of zero skin friction, amethod is needed to integrate the 
boundary layer through the point of emerging reversed flow. 


This method must also account for strong interaction effects 
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ae. to separation and rapid acceleration of the flow 
downstream of the trailing edge. 

The interactive method fulfills these requirements by 
treating the external velocity and displacement thickness as 
unknown quantities. Reflecting the elliptic nature of the 
outer flows, an additional unknown is introduced, but the 
solution can be obtained uSing either the eigenvalue or 
Mechul function methods. 

The Mechul method is preferred as the eigenvalue method 
involves nonlinear problems. In this method the edge 
boundary condition of the direct method is supplemented with 
the interactive boundary condition. The unknown external 
velocity 1s related to its displaced and perturbed 
Senatelons. The unknown functions u(x,y), v(x,y) and ue(x,y) 


are represented in this system of boundary layer equations 





Ou Ov 
a eee (Diev2ecu) 
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where pressure in the y-momentum equation is expressed in 
terms of the external velocity. 

The Mechul function approach assumes’ that’ the external 
velocity, Ue, 1S a function of two arguments, xX and jy, 
allowing for an easy setup of finite difference equations, 
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and avoiding nonlinear eigenvalue techniques. The velocity 
components u and v are required to satisfy the no-slip 
condition at the surface, and u must merge smoothly into the 


edge velocity. 


a& y = OceuiieeO) = OF v(x, Onn = 6 


Uel(X,Ye), Ue(X, Ve) = 


ati Yu ye a VE) 


Uer + — —=t. .— 
Tl ae xXx-& 


1 | d dé 

Uer(X) is the inviscid edge velocity andthe last term, 
the Hilbert integral, approximates the viscosity induced, 
perturbation velocity. 

Interactive methods are useful in both attached and 
separated regions, while direct methods fail at the onset of 
reversed flow, and inverse methods converge poorly. Only at 
the stagnation point singularity are interactive methods 
prohibited. The transformation of the partial differential 
equations into a linear system of algebraic equations is very 
Similar to that of the direct method. The normal coordinate 
y, streamfunction 4, and the external velocity ue are scaled 
with reference to a constant velocity Uo, and the local 


streamwise coordinate x. 


|e (5.25) 
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vx) = ( uo yee £ (oe, 7 ) (o.26) 


Ue (X,Y) 
W(X,n) = ——— (5.27) 
Uo 
Uo is the vector mean velocity. Ue cannot be the reference 
velocity as for Falkner-Skan variables because in this case 
the external velocity is unknown. The first order, semi- 


transformed coordinate system with additional variables U and 


V is 
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and the boundary conditions are 
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The conversion of the flowfield to discrete values is 
very similar to that of the direct method with an orthogonal 
grid, central differences, and two-point averages. In the 
interactive method though, the solution proceeds in the 
downstream direction only. As only downstream disturbances 
are accounted for, backflow causes numerical instabilities. 
A stable integration can be accomplished though, with the 
assumption that backflow velocities are comparatively small. 
The FLARE approximation (Flugge-Lotz and Reyhner), [Ref. 6:p. 
19] sets the streamwise convection term udu/Ox equal to zero 


in regions of backflow. 
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The finite difference equations of the interactive 


boundary layer with the "on-off switch" are now 
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— - Q. (5.34) and (5.35) 


The boundary conditions are also expressed in terms of 


grid or nodal values. A panel method type approximation 
leads to 
2 a 
U; = o, f = 0 
ut = wi wreg +e _ (Won A 
J J’ oJ 1 JJ 5 


where gs and Cas represent a parameter and the diagonal 
element of the interaction matrix, due to a discrete 
approximation to the Hilbert integral. To keep the number of 
generated terms to a minimum, ordinary differential equations 
like the y-momentum equation are centered about the 
downstream face, and partial differential equations like the 
X-momentum equation are centered about the middle of the box. 


The unknowns occur in vectors of four components 


The J quadruplets of unknowns match with 4J equations, 
including 2(J-1) auxiliary relations and (j-1) x-momentum and 


(j-1) y-momentum equations. Each equation corresponds to one 
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of the (j-1) grid rectangles and four boundary conditions. 
The system is linearized around the values of the preceding 


iteration (counter K-1) and a system with Newton iterates 


Sie SUL ON oS, Neo S aera, 
j j j j 
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The overall procedure is a repetitive, linear approach to 
solve the nonlinear system. The numerical solution is again 
obtained using the block elimination method by Keller, except 
that unlike the direct method the vectors of the unknown 


Newton iterates are four dimensional. 
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In matrix-vector form 
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VI. INTERACTION METHODS [Ref. 7:p. 79] 


Interactive methods couple viscous and inviscid flows and 
are intended to compute through regions of flow separation. 
Given their levels of success, these methods have become 
inexpensive alternatives to the Navier-Stokes solvers. 

The simple, classical method of computing the viscous 
flow over an airfoil is: 


ieee tihe velocity distribution is computed for inviscid 
flow. 


2) The inviscid velocity is input to the viscous flow. 


3) The viscous flow is computed by integrating the 
boundary layer equations. 


Now, this method is good at predicting lift and drag, but 
only if the flow remains attached, as information is 
transferred only once from inviscid to viscous regions. For 
more complex flow multiple information transfers are 
required. 

Close coupling is needed to compute flows with separation 
or separation bubbles. A better method than the previously 
outlined classical method for exchanging information between 
viscous and inviscid regions is interaction. The different 
elements of interaction include direct and inverse, inviscid 
and viscous flow solvers. Table 6.1 illustrates the 
different elements. 

The disadvantage of the direct boundary layer method is 
that the equations become Singular at the point of 


5) 


separation. The point of separation may be integrated 
through, however, if the external velocity is computed with a 
predetermined displacement thickness. This method is known 
as the inverse boundary layer method. 

Another problem associated with separation is the 
instability of numerical methods which prohibits downstream 
marching in regions of reversed flow. In the situation where 
flow is reversed the FLARE approximation is used, where the 
momentum transport term udu /Ox is neglected. This 
approximation is not necessarily accurate, but it does allow 
for Continved @calculationcm 

Four interaction models have been developed to calculate 
combined inviscid and viscous’ flows. All procedures solve 
the Laplace equation for inviscid flow and the boundary layer 


equations for viscous flow. The four models are the direct, 


inverse, semi-inverse and viscous-inviscid interaction 
methods. Each model is subject to different boundary 
Conde aonse 


The first three models are considered weak interaction 
methods in that they provide only a loose coupling between 
viscous and inviscid regions. The two regions are treated 
alternately. As indicated in Table 6.1, the viscous flow 
solver calculates the flow in the viscous region and produces 


the boundary condition of the inviscid region. The inviscid 
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TABLE 6.1 INTERACTION ELEMENTS 


BOUNDARY CONDITION 
Flow Inverse 
Poviscid * Zero normal * Prescription of 


velocity at velocity distribution 
the surface 


Viscous No slip * No slip condition 
Congleton 


PReSscriptilon@of * Prescription of dis- 
external velocity placement thickness 
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flow solver calculates the flow in the inviscid region and 
produces the boundary condition of the viscous region. The 
weak interaction methods process either displacement 
thickness or external velocity as an input and =the other 
quantity as an output. 

In contrast, the fourth method, the simultaneous 
interaction method, is considered a strong interaction 
method. A strong method calculates displacement thickness 
and external velocity simultaneously. The foundation of the 


four interaction methods are discussed below. 


A. DIRECT INTERACTION METHOD 

The direct interaction model is composed of direct 
inviscid and viscous flow solvers as indicated in Figure 
6.la. The external velocity distribution is calculated first 
by inviscid computations. The displacement thickness, 65%, is 
then calculated using the external velocity as a boundary 
condition. An updated shape of the displacement body is then 
computed, and all steps are recomputed in order until the 
results converge. As previously stated this method breaks 
down at the point of separation, and is therefore not useable 
in regions of separated flow. However, it is very useful 
where viscous effects are small. The direct method is used 
in the Cebeci Code around the nose and stagnation point of 


the airfoil. 
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Pergtiie. 64.1 Organization of Interaction Methods 
a) Direct, b) Inverse and c) Semi-inverse 
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BB. INVERSE INTERACTION METHOD 

This method was developed to circumvent the singularity 
problems near separation. According to Figure 6.1b it uses 
inverse inviscid and viscous flow solvers. Because of the 
inverse method's very slow convergence, though, it seus 


suitable only at singular points. 


C. SEMI-INVERSE INTERACTION METHOD 

The semi-inverse interaction method incorporates direct 
inviscid and inverse viscous flow solvers such that 
displacement thickness is input to both solvers as shown in 
Figure 6.1c. External velocity is output from both solvers. 
Convergence is ensured with the use of a relaxation formula 


which redefines the displacement thickness distribution. 


On ere oo) = Seameanx) (et + Wy Z 1) ( 6 cil) 


where w 1s a relaxation parameter. 
The numerical weaknesses of the direct and inverse 
methods are improved, but inviscid and viscous regions are 


still loosely coupled. 


D. VISCOUS-INVISCID INTERACTION METHOD 

The viscous-inviscid interaction method ensures a strong 
interaction between the outer, inviscid, and inner, viscous, 
regions. Both the external velocity, ue(xX), and displacement 


thickness, 6*, are unknown quantities. Convergence is 
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ensured through the-interaction law which uses the blowing 
velocity concept. 

The equations are solved through successive sweeps over 
the airfoil surface as indicated in Figure 6.2. For each 
sweep the external velocity for the boundary layer equation 


is written 
Ue(X) = Uer(X) + Sue(xX) (6.2) 
where, 
Uer(X) is the inviscid velocity 
and 


Ue(xX) 1s the perturbation velocity due to the 
boundary layer displacement. 

The perturbation velocity is modeled by the interaction 
law with the help of blowing velocities. The displacement 
effect of a boundary layer is obtained by ejecting fluid at 
the surface of the airfoil as shown in Figure 6.3. 

With a properly arranged blowing velocity source 
distribution on the airfoil surface, the virtual displacement 
body becomes a streamline. 

In determining the source strengths, the displacement- 


body tangential flow condition is represented by 


oe ao~ 





(62.5) 
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Figure 6.2 Viscid/Inviscid Interaction Method 
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Biade surface along which 
sources are distributed 


Pagure: 6.3 The Blowing Velocity Concept 
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Model simplifications are achieved through the use of 
these thin airfoil approximations: 
1) The u-velocity component is considered invariant across 


the boundary layer, as the displacement thickness is 
thin enough to consider the differences negligible. 











2) The blowing velocity, v(x,0), is one half the source 
strength, o(x), on an airfoil represented by a straight 
line. 

O(x) 
= (x,0) 
a 
5” 68v 
=a ( eo*) = — dy 
oO Oy 
dé* Gills 
— ble 6* 
5X dx 
o(x) d 
= —(u.5*) (6.4) 
2 dx 


where (d/dx)(u.6*) is the blowing velocity. 

The blowing velocity once obtained from the source 
strength is then related to the perturbation velocity, 5u.e, 
through the use of the Hilbert integral. 


OU = eee ( Geb) 


2m 





1 | o(&) 


«A pn a 


After substituting equations (6.4) and (6.5) into (6.2) 


the interaction law is obtained. 
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Ue(x) = Uexr(x) + — eo. (6.6) 
Tu Jam GE Ne 


al | d dé 

The numerical implementation of the interaction law 
requires some discrete approximation of the thin airfoil 
integral, equation (6.6). Similar to the panel method, a 


piecewise approximation of the continuous’ blowing velocity 


d(u.65*)/dx allows for piecewise analytical integration. 
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VII. AIRFOIL STUDIES 


Cebeci's interactive aircode was applied to four airfoils 
over a wide range of Reynolds numbers and angles of attack. 
The computer program results were then compared to reported 
experimental data. Unless otherwise stated, 20 iterations 
were used for each computer run, and laminar-to-turbulent 
transitions were determined internal to the program. The 
Significance of the number of iterations will be discussed 


later in the section. 


A. NACA 66;3-018 

Computer results of the NACA 663;3-018 airfoil section were 
compared to the test results of Gault [Ref. 8], which were 
performed in the NASA Ames Research Center 7-by-10-foot wind 
tunnel. The laminated pine model with a 1/8 inch-thick 
mahogany plywood veneer spanned the 7-foot dimension to 
simulate two-dimensional flow. 

Total-and static-pressure surveys, hot-wlire-anemometer 
observations, and detailed pressure-distribution and liquid- 
film measurements were made in regions of separated flow. 
The measurements were obtained for a wide range of angles of 
attack and for Reynolds numbers from 1.5 to 10 million. A 
main purpose of these measurements was to identify locations 


of separation, transition and reattachment. 
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Using the Cebeci Code the 663-018 airfoil shown in Figure 
7.1 was initially tested for section lift coefficients. 
Comparisons were made with Abbott and Doenhoff [Ref. 9] and 
the results are presented in Figures 7.2 and 7.3 for Reynolds 
numbers of 3 and 6 million, respectively. 

Upper surface, laminar to turbulent, transition locations 
are shown in Figures 7.4 and 7.5 for increasing angles of 
attack and for Reynolds numbers of 3 and 6 million, 
respectively. Gault's locations were obtained from pressure 
and hot-wire measurements, which provided near identical 
results. The program transition locations were computed to 
be the point of laminar separation. Note that the transition 
locations shift forward as the angle of attack is increased, 
and they approach the leading edge above 6 degrees. Unless 
otherwise stated, all computer runs7~ used a transition 
constant of GY = 1200. 

Midchord upper surface transitions at less than two 
degrees angle of attack and Reynolds numbers of 1.5 to 10 
Million are shown in Figures 7.6, 7.7, 7.8 and 7.9. In all 
cases the computed predictions were forward of Gault's 
because of laminar separation predicted by the Code. 

While Gault found leading edge separation bubbles, the 
Cebeci Code did not predict them at any angle of attack for 
Reynolds numbers of 3 and 6 million. 

The relationships between separation and transition are 


illustrated in Figure 7.10 for the results of Gault and the 
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Figure 7.3 Lift Coefficient, NACA 665-073, Re= 6 Millies 
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LAMINAR TO TURBULENT TRANSITION 
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Figure 7.10 Upper Surface Transition and Separation, 
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Cebeci code. The experimental results show separation prior 
to transition, whereas the computer results predict 
separation after transition. The importance of this 
difference manifests itself in the difference between the 
measured and computed midchord bubbles, pressure 
distributions and velocity profiles shown in Figures 7.11 to 
wr 

Midchord separation bubbles for angles of attack of 0 and 
2 degrees and Reynolds number of 2 million are plotted in 
Figures 7.11 and 7.12. The lines represent contours where 
u/ue = O. The Cebeci Code midchord bubbles are much smaller 
than those found by Gault. 

Full chord pressure distributions for angles of attack of 
zero and two degrees, and Reynolds numbers of three and six 
Million are shown in Figures 7.13 to 7.16. In each case the 
biggest difference between the experimental results and the 
Cebeci code occurred near the midchord separation bubble 
regions. Figure 7.17 shows a leading edge pressure 
distribution for an angle of attack of six degrees anda 
Reynolds number of three million. 

Midchord velocity profiles are shown in Figures 7.18 
through 7.24 for an angle of attack of zero degrees anda 
Reynolds number of two million, and in Figures 7.25 through 
7.31 for two degrees angle of attack and a Reynolds number of 
two million. These velocity profiles clearly show a big 


difference in bubble sizing. 
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Figure 7.15 Upper Surface Pressure Distribution, NACA 
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Figure 7.16 Upper Surface Pressure Distribution, NACA 
Goa-Cle, AOR = 2°, R = 6 Million 
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Figure 7.19 Upper Surface Velocity Profile, NACA 663-018 
X/C = .62, AOA = 0°, R = 2 Million 
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Figure 7.20 Upper Surface Velocity Profil 
X/C = .64, AOA = 0°, R = 2 Million 
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Figure 7.21 Upper Surface Velocity Profile, NACA 663-018, 
X¥/C = .66, AOA = 0°, R = 2 Million 
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Figure 7.23 Upper Surface Velocity Profile, NACA 66q=@26, 
xX/C =. .69,. AOA = O09 eR = 23M 7on 
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Figure 7.24 Upper Surface Velocity Profile, NACA 663-018, 
Ve = #70, BOA =O. R = 2 Million 
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Figure 7.25 Upper Surface Velocity Profile, NACA 663-018, 


X/C = .58, ACA =52°7 Ree 2 vMaetiron 
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Figure 7.26 Upper Surface Velocity Profile, NACA 663-018, 
mee = ,60  PAOR =e, Ra= 2 Million 
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Figure 7.27 Upper Surface Velocity Profile, NACA 663-018, 
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Figure 7.28 Upper Surface Velocity Profile, NACA 663-018, 


WES t= 64; PROR. =seo, (Re= 2 Mirlion 
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Figure 7.29 Upper Surface Velocity Profile, NACA 663-018, 
X/C°= .665°R0A =92°, Ro] 92 Ma ieron 
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Figure 7.30 Upper Surface Velocity Profile, NACA 662-018 
My) Ge-8.66, ACR = 22) Rk = 2 Million 
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Figure 7.31 Upper Surface Velocity Profile, NACA 663-018, 
By C= 169° AOA = 2° 7 ke  2eMeeeeon 
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In an effort to increase the region of separated flow, 
both the GY transition constant and RHeMmIOCAt ion of upper 
surface transition were adjusted. Table 7.1 shows 
theresults. With a constant of 1200 the transition location 
could not be moved aft. However, when the constant was 
lowered to 200 and below, the transition location, x/c = .69, 
found by Gault, could be used moving the transition inside 
the bubble. While a lowered GY... constant and increased x/c 
transition improved the bubble size, the separation length, 
x/c = .60 to .70, could not quite be met. The best result, 
bubble = .6391 - .7098 (x/c), was obtained with a GY... of 40, 
and an upper surface transition location, XTRU, input of .69 
(x/c). Whether 40 is a suitable value for other foils has 
not been determined. 

Twenty iterations were used on all computer runs for this 
airfoil. To make sure that 20 iterations were sufficient for 
accurate results, Figure 7.32 was obtained. ane. Litt 
coefficient was plotted for each iteration, and as can be 
seen the even iterations produced very minimal changes past 
12. Even between 19 and 20 the change in lift was only 5 x 
10+e.. Therefore, 20 iterations were considered sufficient 


for all computer runs. 


B. NACA 0010 (MODIFIED) 

Similar to the NACA 663-018, computer results of the NACA 
0010 (Modified) airfoil section were compared to the test 
results of Gault [Ref. 8]. The tests were also conducted in 


OW 


TABLE 7.1 EFFECT OF GYTR AND XTRU ON THE LENGTH OF THE 


SEPARATION BUBBLE 


GYTR XTRU (X/7e) SEPARATION (X/C) 


1200 .620792 (1) .6572 - .6929 
5) 2 (2) 
.69 (e2Za) 
400 .69 (2) 
300 .69 (2) 
200 .69 .6391 - .7596 
120 .639164 (1) No Separation 
65 16572 435675 1 
.69 6391 - .7434 
80 . 69 .6391 - .7268 
60 .69 .6391 - .7268 
40 .69 .6391 - .7098 
30 .69 .6572 - .7098 
20 .69 s65725— 27098 
10 . 69 .6572 - .7098 
NACA 663-018 
Reynolds Number = 2,000,000 
Angle of Attack = O Degrees 
Region of Separation, Experimental (Gault) = .60 - .70 (x/c) 


GYTR = Empirical Transition Constant 


XTRU = Upper Surface Begin of Transition, Input or Computer 
Derived 


Lower Surface Begin of Transition = Computer Derived, Each 
Case 


1) Computer Derived 


2) Breakdown in Simulation 
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Figure 7.32 Lift Coefficient Versus Iterations, NACA 
663-018, AOA = 0, R = 2 Million, Transition 
Constant = 1200 


a4 


the 7-by-10-foot wind tunnel at NASA Ames and two-dimensional 
Flow was simulated. 

The 0010 (modified) airfoil shown in Figure 7.33, unlike 
the 663-018, was not computer tested for section 1l1ift 
coefficients as Abbott and Doenhoff [Ref. 9] had no 
comparable airfoil. 

Leading edge, upper surface, laminar to turbulent, 
transition locations were observed, though, and the results 
are shown in Figures 7.34 and 7.35 for Reynolds numbers of 
two and six million, respectively. The Cebeci Code curves 
represent the beginning of transition. 

Full chord pressure distributions for angles of attack of 
zero and three degrees and Reynolds numbers of three and 
eight million are plotted in Figures 7.36 through 7.39. 

Leading edge pressure distributions for angles of attack 
of four, eight and twelve degrees, and Reynolds numbers of 
two and six million are shown in Figures 7.40 through 7.45. 
Of particular interest are the "lump" disparities in Figures 
7.42 through 7.45. A possible explanation for the computer 
program deletions of the lumps is a failure to predict 


leading edge bubbles. 


C. NACA 4412 

Computer results of the NACA 4412 airfoil section were 
compared to the test results of Hastings and Williams [Ref. 
10], which were performed in the 13-by 9-foot low speed wind 
tunnel of the Royal Aircraft Establishment at Bedford. The 
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Figure 7.33 NACA 0010 (Modified) 
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Figure 7.34 Upper Surface Laminar to Turbulent 
Transiticn, Leading Edge, NACA 0010 
(Modified), R = 2 Million 
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Figure 7.36 Upper Surface Pressure Distribution, NACA 
OO10 (Modirited) ,2AGA =56°, R = 3 Million 
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Figure 7.37 Upper Surface Pressure Distribution, NACA 
Q01O0 (Meditied), Bez = 0°. R = 8 Million 
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Figure 7.38 Upper Surface Pressure Distribution, ae. 


0010 (Modified), AOA = 3°, R = 3 Million 
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“igure 7.39 Upper Surface Pressure Distribution, NACA 
0010 (Modified), AOA = 3°, R = 8 Million 
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Figure 7.40 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified), AOA = 4°, 


R = 2° Miron 
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Figure 7.41 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified), AOA = 4°, 
R= 65 Mayon 
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Figure 7.42 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified), AOA = 8°, 
R = 2 Million 
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Figure 7.43 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified), AOA = 89, 
R =869Milli on 
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Figure 7.44 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified), 
AOA = 12°, R= -cereeeron 
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Figure 7.45 Leading Edge Upper Surface Pressure 
Distribution, NACA 0010 (Modified), 
AOA = 12°, R = 6 Million 


sles: 


one meter chord model spanned the full width, 13 feet, to 
Simulate two-dimensional flow. 

Mounted at its quarter-chord point, the model was 
extensively instrumented with static pressure orifices. 
Boundary layer and wake measurements were made at mid-span 
where the 88 pressure orifices were located. 

The main emphasis in the experiment was on defining the 
upper surface boundary layer through separation and into the 
wake. Laser anemometry was used to measure the average 
velocities and Reynolds stresses were measured by hot-wire 
anemometry. 

The 4412 shown in Figure 7.46 was initially computer 
tested with the Cebeci code for momentum thickness. In 
Figure 7.47 the upper and lower surface laminar to turbulent 
transitions were computer derived, xX/Cer upper and lower 
surfaces = 0.00625. In Figure 7.48 the upper and lower 
surface transitions were input as xX/Cer upper surface = 0.01, 
and x/Cer lower surface = 0.11. These values were as close 
as could be input to 0.014 and 0.110 respectively, for the 
downstream ends of the transition trips used in the 
experiment. The differences in transition locations seemed 
to make no difference in computer results. The momentum 
thicknesses still did not agree very well with Hastings and 
Williams' experimental results. 

Figure 7.49 compares lift coefficients from the Cebeci 


code, Hastings and Williams, and Abbott and Doenhoff [Ref. 
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Figure 7.47 Upper Surface Momentum Thickness, NACA 4412, 
AOA = 12.49°, R = 4.17 Million, Computer 
Derived Transitions 
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Figure 7.48 Upper Surface Momentum Thickness, NACA 4412 
AOA = 12.499, R = 4.17 Million, Computer 
Derived Input Transitions 
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10}. The dashed line for the Cebeci code, Reynolds number 
equal to 4.17 million and computer derived transitions, 
Should lie between the curves for Reynolds number equal to 
three and six million from Abbott and Doenhoff. However, it 
lies above the six million curve which further indicates an 
insufficient boundary layer development. With the input 
transition the code prediction for Reynolds number equal to 
4.17 million and angle of attack equal to 12.49 degrees, was 
even higher. Of interesting note though, is that the 
Hastings and Williams prediction, Reynolds number equals 4.17 
fereron and angle of attack equals 12.49 degrees, lies below 
the Abbott and Doenhoff values, possibly indicating an error 
on their part. 

Figures 7.50 through 7.56 show the velocity profiles from 
x/ce = .66 to the trailing edge. In all cases the Reynolds 
number was 4.17 million, the angle of attack was 12.49 
degrees, and the upper and lower transitions were .01 and 
.11, respectively. U/U. indicates the fraction of the 
velocity at the boundary layer edge, and yn/Delta is the 
fraction of boundary layer thickness, where Delta, 6, is 
defined as the layer thickness where the velocity is 99% of 
the edge velocity. 

As these figures indicate, as well as Figure 7.57, a 
Cebeci code velocity profile summation, the code does predict 
separation, but not the extent indicated by Hastings and 


Williams. If the lift coefficient curves, Figure 7.58, can 
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Figure 7.50 Upper Surface Velocity Profile, NACA 4412, 
X/C = .66, AOA = 12.49°, R = 4.17 Million 
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Figure 7.51 Upper Surface Velocity Profile, NACA 4412, 
X¥/C = .74, AOA = 12.49°, R = 4.17 Million 
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Figure 7.52 Upper Surface Velocity Profile, NACA 4412, 
X/C = .78, AOA = 12.49°, R = 4.17 Million 
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Figure 7.54 Upper Surface Velocity Profile, NACA 4412, 
Z/C = .92, AOA =] 912-749° WR 4.17 Miler 
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Figure 7.56 Upper Surface Velocity Profile, NACA 4412, 
X/C = .997, AOA = 12.49°, R = 4.17 Millaen 
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be a reference, then it appears that the Cebeci code predicts 
an underdeveloped flow, too little separation, and Hastings 
and Williams show an overdeveloped flow, too much separation, 
for the given conditions. 

To insure that the Cebeci Code was run correctly, certain 
results by Cebeci, Clark, Chang, Halsey and Lee [Ref. 1] were 
attempted to be duplicated. Figure 7.58 compares two curves 
from Figure 14, [Ref. 1], curves labeled interactive theory 
and interactive theory with a modified transition, with a 
curve obtained using the Cebeci Code with 20 iterations and 
computer derived transitions. Interestingly, the feat and 
third curves of Figure 7.58, interactive theory and Cebeci 
Code, respectively, should be the same, but the two clearly 
are not above nine degrees angle of attack. Even more 
interestingly, the Cebeci Code lift curve in Figure 7.59 
after only 10 iterations does match the interactive theory 
curve . 

Figure 7.60 clearly shows the importance of uSing enough 
iterations to obtain a reasonably accurate solution. 

Finally, Figure 7.61 shows a very good match between the 
pressure coefficients for set conditions of Figure 16, Cebeci 


et al [Ref. 1], and the Cebeci Code, 20 iterations. 


D. FX 63-137 

Computer results of the Wortmann FX 63-137 airfoil were 
compared to the test results of Brendel and Mueller [Ref. 
11], which were conducted in the University of Notre Dame 
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Figure 7.60 Lift Coefficient Versus Iterations, 
NACA 4412, AOA = 12°, R = 1.523 Million 
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.61m xX .61m wind tunnel. Two cast epoxy resin airfoil models 
with chords of .305m and spans of .4m were mounted in the 
center of the test section. Pressure was recorded on one 
model with 96 pressure taps connected through two scanivalves 
to an electronic manometer. Boundary layer velocity 
measurements were obtained on the other model using a 
constant temperature anemometer with a five pm diameter, 
Single-sensor, hot-wire, boundary layer probe. 

Using the Cebeci Code the FX 63-137 airfoil shown in 
Figure 7262 was initially tested horm section. lift 
coefficients with a transition constant of 1200. Reynolds 
numbers of .28, .5 and .7 million were used, and the results, 
shown in Figures 7.63 and 7.64, were compared to those of 
Althaus and Wortmann [Ref. 12]. Interestingly, the Cebeci 
Code predicted low values for Reynolds numbers of .28 and .5 
million, but for .7 million the lift coefficients were nearly 
identical to Althaus and Wortmann up to an angle of attack of 
10 degrees. 

As the purpose of Brendel and Mueller was’ to make 
boundary layer measurements at low Reynolds numbers, a 
computer comparison was unsuccessfully attempted for steady 
flow at a Reynolds number of 100,000 and an angle of attack 
of 7 degrees. With 20 iterations the Cebeci Code failed. 

To understand why the code calculations ceased for this 
case, other computer runs were attempted for the same 


Reynolds number and angle of attack, but with fewer 
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Figure 7.62 Wortmann FX 63-137 
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Figure 7.63 Lift Coefficient, FX 63-137, R = .28 Million, 
and .7 Million 
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iterations. Figures 7.65 and 7.66 show the upper surface 
displacement and momentum thicknesses for steady flow and 
iterations of 2, 4, 6, 8 and 10. As can be seen in both 
figures flow calculations matched very well with experimental 
data up to approximately x/cCc = .55. After that point stall 


occurred and calculations ceased with more than 10 


iterations. Brendel and Mueller experimentally derived 
separation to begin at x/C = .34, but reattachment was shown 
to occur at x/c = .60. Unfortunately the Cebeci Code could 


not predict reattachment for the prescribed conditions. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 


Cebeci's viscous/inviscid interaction program was applied 
to the analysis of steady, two dimensional, incompressible 
flow past four airfoils, the NACA 663-018, 0010 (modified), 
4412 and the Wortmann FX 63-137. Detailed comparisons with 
the available experimental results show that for attached 
flows the essential features are correctly modelled, but that 
Significant discrepancies are found in regions of flow 
separation. These discrepancies are possibly caused by the 
empirical transition modelling used in the present code. 
Future efforts therefore should be directed to the 
incorporation of transition calculations which permit the 
prediction of transition within a separation bubble, such as 
the application of the e”-method proposed by Cebeci [Ref. 
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